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Abstract 

We present three novel forms of the Monge-Ampere equation, which is used, e.g., in image processing and in 
reconstruction of mass transportation in the primordial Universe. The central role in this paper is played by 
our Fourier integral form, for which we establish positivity and sharp bound properties of the kernels. This 
is the basis for the development of a new method for solving numerically the space-periodic Monge-Ampere 
problem in an odd-dimensional space. Convergence is illustrated for a test problem of cosmological type, in 
which a Gaussian distribution of matter is assumed in each localised object, and the right-hand side of the 
Monge-Ampere equation is a sum of such distributions. 
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1. Introduction 

The Monge-Ampere equation (MAE) 

det = /(x) (1) 

is encountered in many areas of numerical analysis and physics, ranging from image processing [H-Q to 
cosmology [JlHl- Here the subscripts Xi denote derivatives in the respective spatial variables; llu^j.^jj^ || is the 
Hessian of u, i.e., the matrix comprised of second derivatives of u. Existence and regularity of its solutions 
was considered in Various strategies were proposed for its numerical solution. A provably convergent 

method for solving the Dirichlet problem for the two-dimensional MAE in bounded convex domains was 
presented in flo|. It is directly linked to the geometric interpretation of the equation that had given an 
opportunity to demonstrate existence of weak solutions [13j |. However, the actual numerical examples in [1Q|^ 
have just up to 25 grid points, which is clearly insufficient to establish the practical feasibility of the methodic 
The MAE can be recast as a minimisation problem; application of algorithms for saddle-point optimisation 
to a two-dimensional MAE was considered in 3) l3| ^ and least-squares minimisation was advocated in [1^ . 



Efficient methods for solution of the discrete optimal transportation problem were developed for application 
to cosmological problems 0-01 ■ It may be perceived that discrete methods correspond well to the physics 
of mass transportation in the Universe — after all, galaxies, as observed by astronomers, are clearly well- 
localised, discrete objects! This argument, however, becomes less persuading, when one recalls that visible 
matter constitutes only a very small fraction of different, "dark" kinds of matter, whose density distribution 
is supposed to be continuous (as opposed to a nearly discrete distribution in clusters of the visible matter). 
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In another group of methods the MAE is treated as a generic nonhnear partial differential e quat ion. 
Application of Galerkin's and finite element methods to a regularised MAE was considered in IT], [iSj. A 
pseudospectral Newton's algorithm was reported to perform well for a MAE in with a smooth r.h.s. (l9l |. 
However, when we implemented this algorithm for a three-dimensional MAE, we found that it failed to 
converge, unless the r.h.s. was a smooth slowly varying function. This is consistent with the observation To'] 
that Newton's method, applied to equations equivalent to the MAE, is useful to improve an approximate 
solution only if the approximation is accurate enough, and easily fails otherwise. Two methods for numerical 
solution of a Dirichlet problem for an elliptic MAE in a two-dimensional convex region U, are examined in 
[20| . The first one employs a finite-difference discretisation of the equation; it is advocated for application 
to the MAE with a possibly singular solution. The second one is an iterative method for the MAE in the 
form of a fixed point problem 



(here V ^ denotes the inverse Laplacian); it is claimed to perform better, when the solution is regular (i.e. 
belongs to the Sobolev space H^|(r2) ). 



A so-called "inexact" iterative Newton-Krylov solver with preconditioning was applied for two- [2lJ and 



three-dimensional 22| grid generation with the properties of equidistribution and minimum distortion. The 



nature of the problem solved in |2lL |22| required only a modest accuracy of solutions, discrepancies the 
order of 10"'^ — 10"'* being acceptable. Examples of computations of such moderately accurate solutions 
illustrating the scalability of the algorithm with respect to the employed spatial resolution were provided 
ibid. 

In Sections 2-4 we derive three alternative forms of (H}: the "second-order divergence" form, the Fourier 
integral form and the "convolution" form, which, to the best of our knowledge, were never presented in the 
literature before. The second form suggests two related methods for computation of space-periodic solutions 
to (HJ. The methods and results of their test applications are presented in Section 6; the Monge- Ampere 
problem of the cosmological type, which we use to test our algorithms, is presented in Section 5, after a 
short statement of the cosmological reconstruction problem. In the Concluding remarks we briefly discuss 
some open questions. 

We note that the MAE's of various kinds are considered in the literature. For the MAE arising in the 
mass transportation problem, typically the r.h.s. of ([IJ depends on the gradient of the unknown function (this 
form of the equation can be found in the original memoir by Ampere j24|). In the theory of PDE's, usually 
one considers the MAE with a known r.h.s. (see 0, [1^). The MAE of this kind arises in the reconstruction 
of the early Universe (see discussion at the beginning of Section 5). Since this cosmological problem is 
the main original motivation for our work, we restrict our discussion to the case of the "mathematician's" 
MAE. However, a straightforward reformulation of our algorithms makes them applicable for the case of the 
r.h.s. depending on the unknown function. 



2. The "second- order divergence" form of the MAE 

The divergence form of the MAE ([1]), in which the l.h.s. of the equation is represented as a sum of the 
first derivatives of certain quantities, is well known (see, e.g., [l3l)- In this section we derive a representation 
of the l.h.s. of the MAE as a sum of second derivatives, which we call the second-order divergence form of 
the MAE. 

Consider a Fourier integral solution 

I u{u3)e"^''du3 

to the MAE in . Substituting the integral into ([T]) and using the identity ioT N x N matrices 



1 

det||ay|| = —J £ii...ijv£ji...jiv JJ^ 



N 

a,; 



where each of the indices ii, i^r, ji, jjv takes the values 1, and ep-^,,,p^ denotes the unit antisym- 
metric tensor of rank N , we find 



det Wu-r 



(-1) 

m 



N 



E 



-Il...IjV^Jl...JJV 



AT 

n 



u{u})uJi^ujj^e'' duj 




(2) 



(3) 



Here a;^ || denotes a iV x iV matrix comprised of columnar vectors a; ^, a;^ . The first factor in 

the integrand of ([Sj is a quadratic function of u), hinting that the l.h.s. of (H]) can be transformed into a sum 
of second derivatives. Indeed, "reverse engineering" of Q reveals an alternative, "second- order divergence" 
form of ([H) in i?^: 

'm E ^^l---^N^3l---jN \y'X^^Xj^---Uxi^_^x,^_^u\ ^ f. (4) 

*1 J ■ ■ ■ J* JV yjl T ■ ■ ,J N 

Equivalence of ((ij and (Hj) is now easily established directly: Using the standard rule for differentiation 
of the products in the l.h.s. of we render it as a sum of products of derivatives of u. Third-order 
derivatives enter such a product only in pairs of the form Uxi^xi^xj^Ux^^xj^^xi^i and hence all products 
involving third-order derivatives cancel out due to antisymmetry of the tensors eii...ijv. Similarly, all terms 
involving fourth-order derivatives cancel out. Consequently, each term in the l.h.s. of ^ gives rise to a 
single term Si^,,,ij^ej-^,,, j^Uxi^xj^---Uxi^xjj^ i and hence the sum is the determinant of the Hessian of u, as 
required. 

This form has an interesting consequence. Suppose, for the sake of simplicity, that a space-periodic 
solution u £ is sought. Let tp he a. smooth function with a finite support (in this context such Lp are 
called test functions). Multiplying the MAE by (f and twice integrating by parts each term in the sum, we 
obtain by virtue of dU 

E £^l...^N^3l■■■3N j ^ ^x^^x,, ■ ■ ■ Ux^^ _^x, ^ ifx^^x, ^ dx ^ J J if d^. (5) 

As usual in the theory of partial differential equations, a weak solution to ([1]), u, can be defined as a function 
satisfying the integral identity ([S]) for any test function (p. (Other definitions of weak solutions to the MAE are 
natural: generalised solutions obtained by geometric construction s |13| a.nd the so-called viscosity solutions 
0; the two definitions are equivalent [23|.) Commonly (e.g., see (ITL l25l|). only one integration by parts is 
performed in this integral identity. Our form is advantageous in that a lesser regularity of the weak solution 
is required for ([5]) to be well-defined. The following argument illustrates this for TV > 2: For any u from 
the Sobolev space W'^_-j^{T^) the integrals in the l.h.s. of ([S]) are well-defined (because by the Sobolev 
embedding theorem this implies Vu £ L2{n-i){T^) and hence u S Loa{T^) ). By contrast, integrals in the 
similar identity obtained by just a single integration by parts are not well-defined for u G W^_i{T^). (Note 
that a viscosity solution for fully nonlinear second-order equations is only required to be continuous [9|.) 
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3. The Fourier integral form of the MAE; positivity and bounds of kernels 



3.1. Derivation 

In the Monge-Ampere-Kantorovich approach to the cosmological reconstruction problems 0, the 
MAE arises for the potential of the inverse Lagrangian map, in which the function / is a ratio between 
matter density at the current epoch and at the much earlier epoch of matter-radiation decoupling (at the 
earlier epoch, the distribution is very close to uniform); thus, / > 0. 

Here, we consider a problem in which the spatial mean of / is positive. Note that for odd iV, if (/) < 0, 
the change of the unknown function u — )• ~u reverses the sign of / and (/) admits the desirable sign (here 
(•) denotes space averaging: 

(/)^ lim / /(x)dx, (6) 

fi^oo \Bii\ Jb^ 

and \Bii\ is the volume of the ball C of radius R). Suppose that u and / have the same spatial 
periodicity; integration of ^ over a periodicity cell then yields 

/dx = 0, 

which is incompatible with (/) ^ 0. Thus wc assume henceforth 

|2 



2 



(7) 



where u' has the periodicity of / and a zero spatial mean (this being just a normalisation), and c is a 
constant. Substitution of Q into ([ij and integration over a periodicity cell yields 



c 



(/>• (8) 



To derive this, note that Ux^xj — c{Sij +u^.^^), where 6ij is the Kronccker symbol, and hence the derivatives 
of u' in the l.h.s. of ([1]) are present only in det and a linear combination of minors of smaller sizes 

of the Hessian of u' . By the algebraic transformation presented in the previous section each such minor can 
be converted into a sum of second derivatives of products of u' with its second derivatives. Thus the spatial 
mean of the l.h.s. of ([T} over any periodicity cell, and hence the mean defined by ([B]), does not involve u' 
and is equal to . 

In a more general formulation, we seek a solution ([7]) satisfying (|5]) and {u') ~ 0, assuming that / and 
u' can be represented as Fourier integrals: 

\/^u' ^ [ 7j{uj)e''^ ''dLJ, u' = ( u'{uj)e''^-''duj, where 51' (w) = -?y(a;)/|a;p, 

(If the problem, defined by (HJ, ([7]) and ([5]), is considered for space-periodic / and u' , integrals in wave 
vectors in what follows are replaced by the respective Fourier sums.) 
Eq. ^ is equivalent to 

detlK.JI =^ / ... / det^ 



' N-1 \ / 7V-1 



Y[ ^(w") U( ^ ) e''^ ''duj\..dLJ^-^du}, 



n=l 
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where ia denotes a unit vector in the direction of a. Our immediate goal is to derive a similar expression 
for the terms in ^ of lower orders in u' . The term of order m is 



E 



T\=m\ri: i„,j„£CT ^ J n: i„ or j„^o- 



Nl 

il,---,iN ,jl,---,jN |cr| = 

(here the sum X]|(t|=t?j. over all subsets cr C {1, ...,-/V} of cardinality m) 

i-iy 



E 



ml 

l<Pl<...<PN-m<N 



X 



X 



W j exp ^^w" • xj du^... du;' 

/ m-l \ 

Y[ viuj") ) ^ - E e''^""duj\..duj"'-'^duj, (9) 



, ri=l 



m — 1 



ri=l 



where 



5] Af,v..,„_(i\...,i'") (10) 



m. 

l<Pl<...<PiV-m<iV 



is the sum of squares of all minors of rank m, 

^^Pl...pN-m{^ )= ^ ^ Sjl ...imPl ...pjV-m (i )im I 

obtained by crossing out rows of numbers pi < ... < pN-m from the N x m matrix 

M = lli^ i™ll 

comprised of m columnar vectors i^, i™. 

Therefore, the problem (IT|) is reduced after the substitution ([7]) to the system of integral equations, which 
we call the Fourier integral form of the problem ^ and (O : 



m=2 ' 

' m — 1 



X 



n=l 



which is now stated in terms of the Fourier coefficients rj(u)) of the Laplacian of the unknown function u' . 
Equations pT|) are valid for all a; 7^ 0; the respective equation for w = is ([5]). 

If the r.h.s. of ([T]) is zero- mean ((/) = 0), the MAE admits the Fourier integral form pT|) (where rj{uj) 
are now the Fourier coefficients of V^it), where the l.h.s. is reduced to a single term for m = iV in the sum 
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3.2. Bounds for the kernels Am in the Fourier integral form 
In this subsection we establish bounds 



Q<Am{i\...,n<—,, (12) 
ml 



provided all vectors have a unit norm. These bounds will play a crucial role for our numerical algorithm. 

Addition to a column of any linear combination of columns i** for s' < s does not change the value 
of any minor A/pj...p„_^ (i^, i™). Using the Gram-Schmidt orthogonalisation process, we change all in 
■Mm to j'^ such that (i) = i^, (ii) for any s, differs from i* by a linear combination of vectors for 
s' < s and thus Am remains unaltered, {Hi) for any s, is orthogonal to all for s' < s. Hence 



where 6s is the angle between i'* and the subspace spanned by {i* | s' < s}. Consequently, 

m 

Am{i\...,n=Am{i\-.,r)l[sm^ Os- 

s=2 

We denote M'm = •■•J™|| and *M'm the transpose of M'^. The identity ^ 

E Af2^...p„_(j':-J'")-det(*X:„A^:„) (13) 

l<Pl<...<pjV-m<A' 

(which does not require orthogonality of {j^, j™}) can be easily proved using the formula 

m 

det 1 1 ay 1 1 = E ^Ji- -i™ IT ' 

We enlarge the set {j'^, j™} by vectors j'^ for s > m to a complete orthonormal basis in i?^. Let U be an 
orthogonal matrix comprised of the N columnar vectors j^, j^, and £ be the N x m matrix, whose all 
entries are except for £ss = 1 for all 1 < s < m. Then Ai'm = and hence 

detCMmM'm) = deiCS^UUE) det(*f£) = dctZ,„ = 1, 

where Tm is the identity matrix of size m. Consequently, 



to! 

s=2 



2 , 



This demonstrates that are sharp bounds. 

3.3. Solution of the MAE for a weakly fluctuating r.h.s. 

If the fluctuating part of the r.h.s. in ([ij, f — (f), is small relative the mean (/), (fTTj) suggests that 
rj{ijj) ~ f{^) and the nonlinear terms in (jlip are small. Then the system can be solved by iteration: 

N 



/m-l \ / m-1 \ 

X H ?/k(^") ]mi^-J2 du;\..du:^-\ (14) 



n=l 
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Theorem 

1°. Suppose positive constants Cq and Ci satisfy the inequality 



E i^51±^ < C„ (15) 
ml 

m=2 



\fiu:)\du:<Co, (16) 



and rjQ satisfies the inequality 



\^K{i^)- fii^)\duj<Ci (17) 



for K — 0. Then iterates |_?^[ ) are globally bounded: |_?7| ) holds true for all K > 0. 

2°. Under the same conditions, the following inequalities are satisfied for all K > 0; 



/ \r]K+iii^) - rjK{t^)\du: < C2 \riK{(^) - r]K-iii^)\du}, (18) 
Ji?" Jr" 

max \rjK -VK{i^)\ <C2 max {rjn (u:) - rjK (19) 



where it is denoted 



Furthermore, if 



N-l 

C2 = 



m— 1 



< 1, (20) 

then iterations \14^ converge to a solution to fill]) , which is unique in the hall |_?7| ). 

The proof is elementary: Inequalities ([T7 1) ~ (|19l) stem from the bounds p^ . By the contraction mapping 
principle, inequalities (IT51) - (pn|) imply convergence of iterations in the norms of Li{R^) and C{R^) in 
the space of Fourier coefficients, yielding a solution satisfying (ITT)) . 

Evidently, equations (|15p and (|20p. where equality is assumed instead of the inequalities, have positive 
solutions for any N. (For instance, Co = -v/S — 4/3, Ci = 1/3 for = 3.) If for the chosen values of Co 
and Ci (|15p holds true, C2 = 1 and the inequality is strict, then all conditions of the Theorem become 
satisfied for a slightly smaller value of Co- 

Note that in view of (1^01) Co < 1 and hence (ITS)) implies / > everywhere. Consequently, our Theorem 
is mostly of interest as a statement about convergence of iterations (fT4|) . Existence of weak solutions was 
proved by geometric methods in for the MAE with an arbitrary positive r.h.s. in a compact convex 
domain. Although our Theorem is significantly weaker because of the strong restrictions on the r.h.s., we 
have proved it for a non-compact domain — the entire space. 



4. The "convolution" form of the MAE 

Following the same algebraic ideas, the MAE can be partially "integrated". We again use the Fourier 
transform of u': 

u' ^ [ u\u3)e"^-''du3, where u'(u;) = (27r)-^ / u'(x)e-*'^ '' dx. 
Jr« Jr« 

Let 
where 
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arg(^(u;)) = arg(M'(a;))/2, if | arg(u'(u;))| < vr; 
arg(f(a;)) = - arg(f(-u;)), if arg(w (a;)) ^ vr. 

That these conditions can be enforced, is elementary for u) ^ 0; the case a; = is not problematic, because 
C(0) = 0. Then 



is a real-valued function. (For a given u' it is not uniquely defined.) 

Let us render the term of order m in u' in the l.h.s. of ([1]) in the terms of ^(x) employing the expression 
®, (Uni), (HSl) and the "identity" 

(27r)-^ / e'^^^'div = (5(x) 
(as usual understood in the sense of generalised functions): 

'm — 1 \ / m — 1 \ 



1 


■ / An ia;i, ■ 


■ •,1a; 









= (-1) 



[ ... f (a;\...,a;") ( fTM'(w") I expfiV w"-x) 

tt)^)" / ... / An ff('^^)w\...,C('^'")'^") expfi Vw" -x I da;i...da;" 



(27r)-^™ / ... / A„ / Va^)e-''^'-^' d^\...J Ve(x)e-''^" ''"dx" 



expl i ^ u;" • X I du}^... du}^ 



n=l 



= 1 f ... /■ det(*||Ve(xi),...,V^(x")||||V^(x-xi),...,V^(x~x™)||)dxi...dx". 

In particular, the linear term (m = 1) is V^u' = /jjjv V^(x^) • V^(x — x^) dx^. 
Therefore, the problem ^ is equivalent, after the substitution JT]) and to 

/ -7 det(*||VC(xi),...,Ve(x'")|l||Ve(x-xi),...,VC(x-x'")|l)dx\..dx™=4- 



(21) 



We call this integral equation the "convolution" form of the MAE. Since it does not involve second-order 
derivatives, it may prove useful for development of an iterative algorithm for numerical solution of the 
problem defined by ([T]), (O and dS]), which involves suitable transformations of the spatial variable. 

If the r.h.s. of (H]) is zero-mean, the MAE also admits the convolution form (jllL where c = 1, the l.h.s. is 
reduced to a single term for m = iV in the sum X]m=i' u'{<jj) in the definition of ^(w) denotes the 
Fourier transform of u. 
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5. A test problem with a cosmological flavour 



One important area of application of the MAE is the reconstruction of the dynamical history of the 
Universe from present observations of the spatial distribution of masses (galaxies, clusters, including their 
dark-matter components). Let us briefly recall the background. Peebles was the first to propose that, from 
the sole knowledge of the current positions of galaxies (from the Local Group which includes our own galaxy) 
without knowledge of their (proper) velocities, reconstruction of the full dynamical history is a meaningful 



goal [27[ . Indeed, the very strong constraint on the distribution of masses at the epoch of decoupling — which 
has to be almost uniform — makes reconstruction a possibly well-posed two-point boundary problem that 
can then be solved by variational techniques. In fact, using convexity techniques, a theorem of uniqueness 
of reconstruction was proved in using the Euler-Poisson equations, which describes the dynamics of 
matter on sufficiently large scales (of the order of a few million light years). In practice, reconstruction 
on such scales was done so far via the Monge-Ampere-Kantorovich (MAK) method [3| (see also [^01) 
which assumes that the Lagrangian map from initial to current mass locations is the gradient of a convex 
potential. This holds exactly for the Zel'dovich approximation [2^ (in which the equation for the velocity 
of matter reduces to the inviscid Burgers equation with a zero r.h.s., and hence the peculiar velocities 
remain constant along trajectories of point masses) and also for its refinement, the first-order Lagrangian 
perturbation approximation (29| . It is then easy, using mass conservation, to derive a MAE for the potential 
of the inverse Lagrangian map (the latter is the Legendre-Fenchel transform of the potential of the direct 
map). By a theorem of Brenier [30| the MAE becomes a Monge-Kantorovich mass transportation problem 
with quadratic cost which, after discretisation, can be solved by optimisation techniques (see [HI for details). 
The r.h.s. of the MAE is equal to the ratio of densities at the present and initial positions. Thus, if the 
initial density was not uniformly distributed, we would have to solve the MAE with the r.h.s. depending 
on the gradient of the unknown potential. However, at the epoch of decoupling (about 380 thousand years 
after the Big Bang) the density of matter was close to uniform, which implies a significant simplification of 
the MAE to be solved: its r.h.s. is known. 

Here we explore for the first time the possibility of directly solving the three-dimensional MAE without 
discretisation on a toy model with a cosmological flavour. 

We assume that the dimension of space is N = 3, and the r.h.s. of the MAE has the following structure: 

The function describes mass distribution for G "objects", which are "localised", if the value of the 
parameter 6 is small compared to the distance between objects. In the cosmological context galaxies or 
clusters of galaxies can be regarded as such objects; then /^^^ describes the total distribution of mass of 
both visible matter and dark matter, in which the galaxies are embedded. Although observations attest, 
that the distribution of visible matter in galaxies is clearly discontinuous, it is astrophysically sound to 
expect that the distribution of all types of matter is smooth due to the prevailing smoothness of the dark 
matter. We assume that objects have density distributions with a Gaussian shape: 

777(9) 

/'^'(r)=(^(^exp(-|r/.(^)n. (23) 

The g-th object of mass m^^^ > is located at r'^^-'. All m*^^\ cr^^' and r*^^^ are independent of the small 
parameter ^ > 0. 

We will be seeking a solution ([7]) with a space-periodic u' to the Monge- Ampere problem ((IJ, where a 
space-periodic r.h.s. / is the sum of "clones" of ([22]) over all periodicity cells. Without loss of generality we 
can assume that the total mass is normalised: 



G 

/ /(r) dr^Y^ ni^a) ^ i 



9 



thus dHD implies c = 1 in ©. 

It is instructive to consider particular solutions to this Monge-Ampere problem. 

5.1. An exact solution to the MAE for a spherically symmetric mass distribution 

In spherical coordinates centred at r*^^), ([T]) becomes, for spherically symmetric u and /, 

^d^u ( duY _ 3 



where p—\v — r'-^-'|. This equation has an obvious solution 

u{p, 5) = j'^ £ " r^fir) dr^ dr' , (24) 
whose first derivatives are uniformly bounded: 

the second derivatives are 0{S^^): 

2 / p/s \ ^^^^ / \ f p/^ ^ 



This estimate follows from the following inequalities: 

• \XiXm/p'^\ < 1; 

• for p < 6, f{p/S) is uniformly bounded and 



rp/S 

Q.{p/5f < I r^f{r)dr < c{p/5f 







for some positive constant c and c; 

• for p > (5, {p/5Y f{p/5) is uniformly bounded and 



< c' < / r^f{r)dr < c' 

for some constant c' and c'. 

Moreover, if p is larger than a positive constant, f{p/6) = o{S^), and hence all second derivatives of u 
are uniformly bounded outside any sphere of a fixed radius. 

The properties of the solution discussed in this subsection are implied just by a fast decay of / at infinity 
and the spherical symmetry of the r.h.s.; the assumption that the profile is Gaussian is unnecessary. 

5.2. A one-cell solution for G > 1 spherically symmetric localised objects 

A solution for G > 1 objects can be expected to admit a power series expansion 



*(r) = ^^.„(r,<5)<5", 



ri=0 



where (by analogy with (HH)) each M„(r, S) and its first derivatives are 0{S°), and the second derivatives are 
0(<5-i). 
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Naively one could expect interaction of localised objects to be asymptotically unimportant, and hence 
to obtain in the leading order the sum of individual one-object solutions: 

Mo(r) =^u(9)(|r-r(f)|,(5). (25) 

9=0 

Here M^^-'(|r — r^^)|,(5) is the one-object solution ([M)) for the g-th object. Let us inspect, to what extent such 
a conjecture might be true. 

Consider a neighbourhood of an object 7. In the leading order, the l.h.s. of ((T|) is det ||woa;;a;J|. Substi- 
tution of (|25p into the determinant yields a sum of triple products of second derivatives of various m^^^. As 
shown in the previous subsection, for 5^7 each second derivative of 0(5'^) in this neighbourhood 

(the distance between the objects g and 7 is a fixed positive constant). Thus, only triple products of second 
derivatives of w*^^^ contribute terms of the leading order, 0{S^'^), the one of the l.h.s. These products con- 
stitute det and thus match the term in the r.h.s. Nonlinear interaction of pairs of solutions 
for individual objects is at a lower order, 0{S~'^), and that of triplets of one-object solutions - still weaker, 
0{S^^). Hence, after the finite sum is substituted into the MAE, the highest order terms do cancel 
out, and in this respect the naive conjecture is confirmed. 

Nevertheless, the one-cell solution does not represent the leading order term of the solution to the 
problem at hand. To see this, note that at large distances ((25|) exhibits a linear growth in |r|, and not 
a quadratic one, as the form of the solution ([T]) requires. Furthermore, after subtraction of the quadratic 
profile c|rp, the remaining part u' must be space-periodic, which is clearly not the case for (l^5t . We cannot 
use (PSl) to construct a global solution conformant with this periodicity requirement by the procedure, used 
to obtain the periodic r.h.s. / from the individual density profiles ([22]). because the sum of periodically 
distributed "clones" of one-cell solutions (P5|) is infinite. Moreover, for such hypothetical space-periodic sum 
there would be infinitely many pairwise interactions between objects yielding products of the order 0((5~^); 
hence, even if the sum of such products is finite, it will not necessarily remain 0{S~^). 

We observe that the pairwise interaction does not become weaker when the distance between the inter- 
acting objects grows. Thus, even in the high-contrast limit the solution that we are seeking does not reduce 
to individual interactions between objects; their interaction is essentially collective, which makes it hard to 
predict the structure of the solution. 

6. Computation of a solution to the MAE v^rith a space-periodic r.h.s. 

In this section we present two iterative algorithms for numerical solution of the space-periodic Monge- 
Ampere problem. One version (AICDM), employing numerical improvement of convexity and discrepancy 
minimisation stabilising the iterative process, is suitable for computation of solutions for everywhere positive 
right-hand sides / (see Subsection 6.4). Another one (ACPDM), involving continuation in parameter and 
discrepancy minimisation, does not require this condition to be satisfied (see Subsection 6.2). They both 
rely on the basic algorithm for iterative solution of a fixed point problem for the MAE (see Subsection 6.1). 
Test applications of the two algorithms are considered in Subsections 6.3 and 6.4 . We assume (/) 7^ 0, 
however, reformulation of ACPDM for the case of a zero-mean space-periodic r.h.s. is straightforward. 

6.1. A basic solver (FPAR) 

If the amplitude of fiuctuation of / is small compared to the mean, — more precisely, if (jl6l) is satisfied, 
— the Theorem (see Subsection 3.3) establishes existence of a solution to (HI, that is a perturbation of 
c|xp/2. If, in addition, (^0)) holds true, the Theorem guarantees that iterations defined by (ITU) converge 
to the Laplacian of the perturbation; this offers an algorithm for numerical solution of the MAE. If either 
of the conditions ([T51) or (pn|) are violated (which is the practically interesting case), iterations ([T^ do not 
necessarily converge. Different algorithms are necessary, taking this into account. 
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Suppose the kernels Am in (jlip are frozen and take constant values am, respectively. Then the new 
system of equations (jlip can be easily solved, since it is the Fourier form of the polynomial equation 

N 
m=0 

in the physical space. This observation suggests the following algorithm. We express pT|) in the terms of 
rj = V^u' as 

N 

Y,amir^fl{f)+F{T^), (26) 

m— 

where 

N 

F{7j) = 1 + ^ flmr," - det + <5, 

■m— 1 

and implement iterations 



^■y II 1 



N 



J2''rnVK=f/{f)+F{VK-l) (27) 



n=0 



in the physical space. 

The applicability of this algorithm depends crucially on the choice of the coefficients am ■ In view of the 
positivity and the bounds for the kernels (jl2p . we impose 



< am < — r for to > 0. 



Furthermore, it seems practical to set 



ai = 1, am = Tr~l fo'^ TO > 1, (28) 
because for ai = 1 the linear in u' term is treated exactly, 

Fiv)^y2 - / (A„,(i,^i,...,ia;"0-«m) n'?(^")N'^PrI]'^"-^H'^'---^'^'" (29) 
and the median values (PSI) of Am for to > 1 minimise the ranges of the kernels 

in (gni). Since for the coefficients (gS]) 

|^m(iu;i,---,ia;"«) -flml < flm, m>2 

the l.h.s. of (P7|) "captures" the nonlinear behaviour of the l.h.s. of (dJ), the algorithm has chances to 
converge. However, for the extreme values of Am the values of the kernels in ([^^ are as large as the 
respective medians, and therefore convergence of iterations (j27p . (051) is not guaranteed. 

At least two strategies can be proposed for setting the value of ao (which is a free parameter in the sense 
that is not required to be satisfied for u = 0): 

1°. At each iteration oq is tuned, so that (77) — (which holds true for rj — V^u'). 
2°. We set oq = 0. 

Note that although in the subspace of space-periodic functions the inverse Laplacian is defined for zero- 
mean scalar fields only, (rjn) is not required to vanish, because the mean is removed by differentiation in 
the Hessian before the inverse Laplacian is evaluated. 
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For any odd N, the equation (|27p in r/K with the coefficients ([25]) has a unique root for any r.h.s. To 
check this, it is enough to establish that the derivative 0^(7]) of the l.h.s. of (P7|) is positive for any ry. For 
the choice of coefficients , 

D^M - D'j^ir^) - ^1 + jfzjy) ■ (30) 

At a minimum D'j^{r]) = and hence (15(11) implies that at the minimum Djv > 0, proving monotonicity of 
the l.h.s. of (|27l) . It can be shown similarly, that for any even N the number of roots is or 2. Thus, the 
algorithm is guaranteed to be applicable for odd N only. 

It can be proved, that the Theorem also applies for iterations (1271) (with the same constant Ci bounding 
the same norm of solutions). 

In an application to a test problem inspired by cosmology (see Subsection 6.3), this algorithm, which 
we call FPAR (Fixed Point Algorithm for the Regular part of the MAE) produces a sequence of iterations, 
initially converging, but subsequently blowing up: linear instability sets in, and hence the respective unstable 
mode must be removed. 

6.2. A more advanced solver (ACPDM) 

The behaviour of the basic algorithm FPAR suggests that it should be embedded as an engine within a 
more advanced algorithm. In this subsection we present such an advanced solver, ACPDM [Algorithm with 
Continuation in a Parameter and Discrepancy Minimisation) . 

Consider a generalisation of (^5]) : 

Q{W,P) = Q, (31) 

where 

N 

Q{rj;p) = J2 «™'7'" ' f/if) ~ pHv)- 

m=0 

Here the coefficients (j28p are assumed, and the new parameter p is confined to the interval [0, 1]. For p = 0, 
([5T|) is just a set of polynomial equations of degree N\ for p = 1 it reduces to (|^5|) which is equivalent to 
([T]). Continuation in the parameter p is implemented: (1311) is solved for a set of values Pj, increasing from 
to 1, and an initial approximation of the solution for p = pj is obtained by polynomial extrapolation of 
solutions for all pj' < pj. (Numerical extrapolation requires performing quadruple precision computations, 
if the number of nodes exceeds roughly a dozen.) For any p, solution of pip involves iterations similar to 

N 

Y^a^i^^^ f/{f)+pF{r^K-i). (32) 

m=0 

The r.m.s. discrepancy 

d{ri) ^ ^{{Q{v;p)-{Q{mp))Y) 

is used in the termination condition. Note that {Q{ri; 1)) — and hence 

div) ^ .\/((det||V-2r?,,,^.+5.,|| -//(/) )2) 

for p = 1; in particular, d{rj) — for p = 1, if and only if for u' = V^^ry and the normalisation ([5]) the field 
dZl) is a solution to the MAE. 

Let (•, •) denote a scalar product and || • || the induced norm of a scalar field: ||v|| = y^v/v)- In. the test 
runs reported in the next subsection, the scalar product of the functional Lebesgue space L'^{T^) has been 
assumed: 

(u, v) ~ / w(x)w(x)(ix. (33) 

Other products (with different weight functions introduced in the above integral) have been also considered 
(see Subsection 6.4). 
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If 7]' and 77" are approximate solutions to the generalised MAE (|3ip . then 



Q{v';p) - Q{v";p) = A{fj' - rf) + 0{\W - v''^), (34) 



where A is the linearisation of pip around the solution. Consequently, the concept of minimisation of the 
residual in Krylov spaces and approaches for its realisation can be borrowed from solvers for linear prob- 
lems (such as the Generalised Conjugate Gradients Method |3l|). ACPDM involves sequences of stabilised 
iterations described below, which exploit this concept. To make such an iteration, the following data is 
required: an approximate solution r]K, and two sets of S scalar fields, Vs{k) and Ws(x), < s < S, where 
< 5* < 5'niax (as well as some other quantities computed at the previous stabilised iteration), and S'max is 
a parameter of the algorithm. It is supposed that all Vg are mutually orthogonal with respect to the scalar 
product (•,•), and 

V, ^ Aws + 0{\\wsf) (35) 

for small Wg- 

A sequence of stabilised iterations of A CPDM is initialised using the current approximation rjo by com- 
puting rji — rj'^ as a solution to ([5^ ior K = 1, and setting S — (i.e. the sets Ws{x) and Ws(x) are empty). 
For K > 1, a stabilised iteration of ACPDM consists of the following steps: 

i. At each grid point in the physical space solve equation ([5^ : 

N 
m=0 

ii. Compute Q(?7^). 

Hi. Orthogonalise v' = Q{rj'j^) — Q{'^'k-i) all 1 < s < S , with respect to the scalar product (•, •), and 



set 



iv. Compute 



s—1 ^ ' s—1 ^ ' 



5+1 



{Q{v'k).^s), 



V. If 5 < S'max, increase S by 1; otherwise (i.e., if S* = S'max), discard vi and wi, and decrease by 1 the 
indices s of the remaining Smax fields Vs and Wg ■ 

vi. Compute the r.h.s. of p2p for r]K+i substituted in place of r]K-i, the field Q{r]K+i) and d{rjK+i). If 
d{rjK+i) is less than a given small threshold, then r]K+i is the desired approximate solution and computation 
for the present p is finished. If ||(5(r7ii:+i)|| > ||Q(»/a:)||, then the current sequence is terminated. 

A few comments are in order. Clearly, vs+i and ws+i obtained in step Hi, possess the required properties: 
vs+i is (•,•)— orthogonal to all Vs for s < S, and ((35|) holds for vs+i and ws+i by virtue of (p4| . The 
coefficients in the sum computed in step iv minimise the discrepancy ||Q(?7^ — X]s=i 9sWs)|j, assuming 
all QsWs are small and hence their quadratic contributions are negligible. The minimisation plays a dual 
role: on the one hand, it stabilises basic iterations (j32p . removing the instability modes as soon as they 
become substantial; on the other, it significantly increases the efficiency of the algorithm (up to a factor 20 
compared to other algorithms relying on iterations p2l) ). The assumption that nonlinear terms are small 
can be incorrect; also, as a result of accumulation of the neglected nonlinear terms after a number of steps, 
at some stage Vs can cease to approximate Aws accurately enough. If the inequality ||Q(?/x+i)|| > ||Q(^a:)| 
is found to hold true in step wi, we interpret this as an indication that the adverse effect of nonlinearity 
has become significant, and then the algorithm breaks the current sequence. In order to reduce the adverse 
influence of nonlinearity, a small number S'max may be chosen; in our runs, S'max = 5. 

Now we can assemble the algorithm from the building blocks, discussed above. ACPDM performs 
continuation in the parameter p. For a given p, it starts by carrying out basic iterations p2p . As soon 
as convergence slows down (we have used the condition d'^{riK) > d'^ijjK-i) /'^), the algorithm switches 
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to perform a sequence of stabilised iterations. If the sequence terminates in step vi because nonlinear 
effects became significant, the algorithm proceeds by performing basic iterations (|32p . Usually it takes a 
small number of them for the instability to set in, and as soon as ACPDM detects that the inequality 
d^iVK) > ot<P{riK-i) holds true, it starts a new sequence. Here a is a parameter, which can slightly exceed 
1 (in order to allow the transients to die off and the dominant instability modes to set in, so that the latter 
could be efficiently removed in subsequent stabilised iterations); in our test runs we have used a = 1.05 . 

6.3. Application to a test problem inspired by cosmology 

We have chosen to employ a problem of the kind discussed in Section 5 as a test-bed for our algorithm, 
because cosmological applications of the MAE are probably the most important ones. The positiveness of 
the r.h.s. (implying convexity of solutions [13;]) is not required for application of the algorithm. 

We are seeking a solution ([7]) with a space-periodic u' to the Monge-Ampere problem ([ij , where a space- 
periodic r.h.s. / is a sum of "clones" of ([2^ over all periodicity cells. Our poor man's Universe involves 
G = 3 objects in the periodicity cell T^{[0, 1]^), centred at the origin: —1/2 < Xi < 1/2. They are described 
by the following parameter values in psp: 



■ ^(2) . ^(3) ^1 1 1 



s = i, 



18 ' 27 ' 32' 



r^'^- r(2^ = i(l,-l,-l), r(3' = i(-l,-l,l). (36) 

Noting that /(x) achieves its maximum at r*^'^^ and its minimum at -(1, 1, 1), we can estimate the contrast 
number of the problem as 

— ^ ^ w 1.313 X 10^ 

4(12cxp(-18) + 8exp(-18) + 16exp(-32)) 

(the factor 4 in the denominator accounts for the replicas of , located in the neighbour periodicity boxes 
at a distance 1 / -^72) . 

Iterations defined by (j32]) have been performed on a uniform grid comprised of 64^ points in the periodicity 
cell T'^([0, 1]"^). We have evaluated det +Sij \\ by the pseudospectral method with dealiasing (for N ~ 3 

this requires computation of the derivatives u'^.^j the physical space on the twice finer grid involving 128'^ 
points). Apparently dealiasing does not play any important role in these computations; this is consistent 
with the fast decay of the energy spectrum of 77 = V^w' by 11 orders of magnitude. 

We have tested both strategies for choosing gq (see Subsection 6.1), as well as a hybrid strategy, in which 
ao ~ 0, but after solving (j27p we reset rjK '.— rjx — (vk)- It turns out that the hybrid strategy and 1° 
are very close in the number of iterations necessary to obtain a solution to the MAE to the same accuracy 
d{riK) < 10~^°. However, this implies that the hybrid strategy is several times more efficient in the terms of 
CPU time, since it requires just one evaluation of solutions to (j27p . while evaluation of oq following strategy 
1° requires several such evaluations. The strategy 2° has proved to be the most efficient one, both in the 
terms of CPU time and the number of required iterations. (The same nodes pj were used in all the test 
runs.) In the remaining part of the subsection we discuss convergence of the advanced method with ao = 0. 

ACPDM performs a combination of basic iterations (|32|) and stabilised iterations, which have different 
computational costs. The most time-consuming operation, determining the cost, is computation of the 
determinant of the Hessian in F(ri). There are two such operations in a stabilised iteration, and a single 
one in a basic iteration. Consequently, the former is approximately twice longer than the latter. To enable 
an accurate comparison of performance of various versions of our code, we report the number of evaluations 
of determinant in each run (a run is terminated when the obtained approximate solution 77 satisfies the 
accuracy requirement ^(77) = 10^^''). Note however, that the main bulk of computations is performed by 
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Figure 1: Number of evaluations of tlie determinant of tlie Hessian (vertical axis, logarithmic scale) performed by ACPDM in 
successive computations of approximate solutions r;(pj) to the generalised MAE II31II . satisfying d(»?(pj )) < 10"^", for numerical 
solution of the test MAE I I31I1 (see Subsection 6.3) Horizontal axis: the index j numbering consecutive nodes pj in the mesh 
II37II . The initial approximation for a pj is obtained by the polynomial extrapolation of solutions for pj/ with j' < j. 



doing stabilised iterations; the number of basic iterations is usually below 1% of the number of stabilised 
ones. 

In all runs, the initial iteration is a solution to (|5T|) for p = 0. If ACPDM is applied to this field for 
p = 1, iterations quickly become chaotic and cease to converge (this has necessitated to include into the 
algorithm continuation in the parameter p). Initially, ACPDM has been applied for uniformly distributed 
nodes pj — j / J, j ~ 0, J— 1 for J = 20, for which the algorithm has shown a remarkably fast convergence 
(see Fig. 1). A polynomial 20- node extrapolation yields an approximate solution to our test problem to the 
accuracy d{r]) = 0.47 X 10^2 ^^(^^ = o.02, where 

dooiv) = max det + 5ij\\ - J/(f) 

When it is used as an initial approximation for a run for p = 1, convergence is slow and the pattern of 
convergence is erratic. It takes 38 685 evaluations of the determinant of the Hessian for the two discrepancy 
norms to decrease to 0.99 x 10~^ and 0.39 x 10~^, respectively. At this stage convergence stalls. We have 
therefore got into a local minimum of d{r]), out of which no exit can be found; we will refer to it as a spurious 
minimum solution. 

A better approximation to the solution ^ to our test problem is obtained by adding more nodes near 
the right endpoint p = 1. We have chosen to add J' ~ 13 nodes, the complete p-mesh being 

Pj=j/J, .? = 0,...,J-1; pj+^_, = l-i2^J)-\ J = !,...,.]'■ Pj+r = l (37) 

(information on convergence at the new nodes Pj+j is also included in Fig. 1). Initial approximations 
at each new Pj+j become more and more accurate in the interval < p < 0.55 {j = 11), then the 
discrepancy d{rf) of the initial approximations starts growing and admits the maximum 0.52 x 10"'^ for 
p = 0.9875 (j = 21), and subsequently decreases again. The rate of convergence progressively falls down as 
Pjj^j approaches 1: the number of evaluations of the determinant of the Hessian yielding solutions of the 
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Table 1: Number of evaluations (NE) of the determinant of the Hessian, performed by AICDM in the course of computation 
of approximate solutions Q of varying accuracy to the test MAE with the r.h.s. defined by II22II . II23II and II36II . discrepancy 
dociv) ior these solutions, and wallclock duration of runs (DR, seconds) for computation of approximate solutions to the test 
MAE by E AICDM. 





10-3 


10-4 


10-5 


10-*^ 


10-^ 


10-8 


10-9 


10-10 


dociv) 


0.64x10-2 


0.86x10-3 


1.02x10-3 


0.23 xlO-'' 


0.28x10-5 


0.39x10-*^ 


0.36 xlO-'^ 


0.41x10-8 


NE 


64 


126 


251 


439 


1018 


1515 


2 030 


2 653 


DR 


3.38 


4.69 


11.17 


42.7 


81 


203 


435 


774 



Table 2: Number of evaluations (NE) of the determinant of the Hessian, performed by AICDM with the use of various 
weight functions II38I I in the scalar product (■, ■). All runs are terminated as soon as the accuracy d{r]) < 10"^'' is obtained; 
d(jj) = 0.66 X IQ-^" for q = -1, d{'q) = 0.95 X 10-^° for q = -3/4 and d{r]) = 1.0 X 10-^'' in all remaining cases. 



q 


-1 


-3/4 


-1/2 


-1/4 





1/4 


NE 


3169 


2 619 


2 465 


3 743 


2 653 


2 571 


dooiv) 


0.09 X 10-8 


0.14 X 10-8 


0.21 X 10-8 


0.62 X 10-8 


0.41 X 10-8 


0.39 X 10-8 



q 


1/2 


3/4 


1 


5/4 


3/2 


7/4 


2 


NE 


2 568 


2 643 


2 523 


2 489 


2 481 


2 505 


2 526 


dooiv) 


0.37 X 10-8 


0.39 X 10-8 


0.36 X 10-8 


0.38 X 10-8 


0.39 X 10-8 


0.41 X 10-8 


0.37 X 10-8 



desired accuracy d{r]{pj^j)) < 10 markedly increases. However, ACPDM does not stall any more. A 
polynomial extrapolation involving the 33 nodes delivers an approximate solution u' for p = 1 to the accuracy 
d{rj) = 0.78 X 10-"^ and dooiv) = 0.20 x 10-^, and after further 1 885 evaluations of the determinant ACPDM 
yields an approximate solution with the two discrepancy norms down to lO-^" and 0.26 x 10-8, respectively. 
In total, 9 814 evaluations of the determinant of the Hessian are involved in computations on the mesh p7l) 
providing the solution to the MAE. 

The geometry of the obtained solution — isosurfaces of u' and V^w' — is shown in Figs. 2 and 3, 
respectively. The structures disclosed by the isosurfaces of V^m' (Fig. 3) are clearly associated with the three 
objects incorporated into the r.h.s. ([^^ of the test problem. The figures also reveal a subtle interaction of 
the objects along the lines connecting their centres (see in Fig. 3 the tube-like structures, which connect the 
regions of higher values of V^w' encompassing the centres of objects). 

Figures 4 and 5 display the structure of the field Vm', determining the displacement of mass in the test 
Universe. It turns out that u' possesses hidden symmetries, namely, mirror reflection symmetries about any 
plane that is parallel to a coordinate plane and contains a pair of objects. (Deliberately, our computations 
are not sped up by exploiting these symmetries.) Because of the symmetries, the component of Vit' normal 
to such a plane vanishes on this plane. To illustrate the behaviour of the components of the gradient not 
shown in Fig. 4, we present in Fig. 5 their isolines on coordinate planes, which pass through the centre of the 
periodicity cube T3([0, 1]3) and hence are displaced from the planes shown in Fig. 4 by a quarter of period. 

6.4- A solver for MAE with an everywhere positive r.h.s. (AICDM) 

Clearly, it is desirable to avoid the refinement of solutions for nodes 0.95 < pj < 1, which has proved 
necessary in application of ACPDM to our test problem. Comparison of the spurious minimum solution, 
obtained with the 20 nodes pj, and the more accurate one is instructive. The maximum discrepancy between 
the two approximations of r/ = V^m' for these solutions, equal to 0.132, is attained at the minimum of the 
r.h.s., (1/4, 1/4, 1/4), and all points, where the discrepancy is larger than 0.02 are located within the distance 
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(-0.5,-0.5,-0.5) 



Figure 2: Isosurfaces of the solution Q to tiie test MAE (presented in Subsection 6.3) at tlie levels of a half and 1/8 of the 
maximum. The periodicity cell T3([0, 1]^) of u' is shown. 
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Figure 4: Vti' for the solution to the test MAE (presented in Subsection 6.3) on cross sections of the periodicity cell T^([0, 1]"^) 
that are parallel to coordinate planes and contain pairs of objects: X2 = —1/4 (a), xi = —1/4 (b), 0:3 = —1/4 (c). (Due to the 
symmetry of u' about each of the three planes, components of gradients normal to the planes are zero.) The labels Xi refer to 
the Cartesian coordinate axes, parallel to sides of the cross sections. Stars show locations of the three localised objects II36I I on 
the cross sections. Gray-scaling reflects the masses of the objects (black, gray and white stars; the objects at r^^\ g = 1,2,3, 
respectively). 
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Figure 5: Isolines step 0.02 of normal components of Vu' for the solution to the test MAE (dashed lines: negative values, solid 
lines; zero and positive values) on Cartesian coordinate planes X2 = (a), xi = (b), 0:3 = (c). The labels Xi refer to the 
Cartesian coordinate axes, parallel to sides of the shown cross sections of the periodicity cell T^([0, 1]^). 
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1/16 from this point. The Hessian of ([7|) computed for the spurious minimum solution is negative at this 
point, i.e. this sohition ([7]) is not a convex function, as it has to be [l^. This has suggested to develop a 
modification of ACPDM, presented here, which we call AICDM (Algorithm with Improvement of Convexity 
and Discrepancy Minimisation) . 

The algorithm proceeds as an extension of the ACPDM operating with a single node p — 1; computa- 
tion of a good initial approximation by extrapolation in p is unnecessary. AICDM involves an additional 
procedure: improvement of convexity of an approximate solution, which is performed once the condition 
d'^iVK) < PfPiVK') is satisfied. Here K is the number of the current iteration, K' is the number of the 
iteration, at which the last previous improvement of convexity was performed, and /3 < 1 is a constant factor 
(we have chosen (3 — 0.01 in the test computation reported in this subsection). 

Improvement of convexity in AICDM could be performed by computation of a convex hull of the iterate, 
but we prefer to apply a numerically simpler procedure. At each grid point in the physical space, we compute 
the eigenvalues A of the Hessian of u' = V^^{rjK — {vk)) (they are all real, since the Hessian is a symmetric 
matrix). If they are all larger than —1, then the approximate solution given by Q for the current iterate rjK 
is locally convex at the point under consideration (in this discussion we assume that, after normalisation, 
c = 1). Accordingly, if the minimum eigenvalue Amin exceeds —1, no action is taken. Suppose now the 
contrary, i.e. Amin < — 1- The minimum eigenvalue would become —1, if at this point each second derivative 
wf^.j,. is increased by —1 — Amin, and hence the Laplacian r]K — V^u' is increased by 3(— 1 — Amin)- Following 
this observation, at the points where Amin < ^1 we increase r]K by 6(— 1 — Amin) (by choosing the factor 
6 instead of 3 we are "overimproving" tik)- It is not guaranteed, of course, that |rp/2 + V~^(?7^ — (77^)) 
is convex for the resultant 77^, since this procedure changes all mixed derivatives at each point and does 
not guarantee that each second derivative u'^.,j.. at the points of local non-convexity is increased by the 
same amount, or that non-convexity does not appear at new grid points. Hence we proceed, repeating the 
procedure till each eigenvalue at each grid point becomes larger than —1 — (i(?7x)/2. In our experience, only 
a small number of such iterations is necessary to enforce this condition. The quantities 6(— 1 — Amin) are 
small, being at most comparable with the global discrepancy d{rjK)\ although the discrepancy d{ri'j^) for the 
"improved" approximation 77^ can exceed the discrepancy d{rjK) for the original approximation rjK, it thus 
turns out that the increase (compared to d{rjK)) is usually modest. 

An approximate solution to the test problem (formulated in Subsection 6.3), satisfying ^(77) — 10"^*^ 
and drxi{ri) — 0.41 x lO^'*, has been obtained by the AICDM. This has required 2 653 evaluations of the 
determinant of the Hessian. Table 1 illustrates the deterioration of convergence in this run, as better accuracy 
numerical solutions are successively found. Note that computational cost of each iteration in improvement of 
convexity of an approximate solution slightly exceeds that of computation of the determinant of the Hessian, 
and each iteration is included into evaluation counts (one evaluation per a convexity improving iteration) 
presented in Tables 1 and 2. 

In this run, as in runs reported in Subsection 6.3, the Lebesgue space scalar product p3|) has been 
assumed. We have also inspected convergence of AICDM employing scalar products 



(see Table 2). For q > 0, more prominence is given to discrepancy in the regions, where the r.h.s. of the 
MAE, /, admits relatively high values; for g < 0, the opposite happens, i.e. discrepancy in the regions, 
where the values of / are relatively low, is given more weight. Duration of the shortest run, for q = —1/2, 
with a sequential code on a 3.16 GHz Intel Core Duo processor is 25 min. 19 sec. The computations do 
not reveal any clear dependence of the efficiency of AICDM on the power q — in a large interval of q the 
variation of the number of iterations is just several per cent. The low sensitivity to the value of g > is 
probably linked to the relative smallness of the region where / > (/). 




with weight functions 




(38) 
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6.5. An enhanced version of AICDM 

The efficiency of AICDM can be further improved without employing new mathematical ideas. For 
instance, we have explored the possibility of gradual refinement of approximate solutions with increasing 
spatial resolution. The test problem, considered in the two previous subsections, was solved with the 
resolution of (16M)^ Fourier harmonics, where the integer M is successively increased from 1 to 4. When 
seeking a solution with the resolution 64^ harmonics, which satisfies the accuracy requirement c?(ry) < 10"'^, 
the intermediate-stage computations with the resolution of (16Af)"^ harmonics are terminated as soon as 
AICDM finds a solution to the accuracy 

d{Tj) < io-™M«,2.5M)^ 

We call this algorithm enhanced AICDM, or EAICDM. Durations of runs for integer k > 3 are shown 
in Table 2 (we do not present total numbers of evaluations of the determinant of the Hessian in runs by 
EAICDM, since the times of their computations vary significantly with the resolution). 

A further acceleration is likely to be achieved by avoiding dealiasing, but we did not explore this possi- 
bility. 

7. Concluding remarks 

We have presented new forms of the Monge-Ampere equation in : the second-order divergence (|4]), 
Fourier integral (jlip and convolution (I21[) forms. They have been derived under the assumption that the 
MAE ([T]) has a r.h.s. with a non- vanishing spatial mean, and hence ([T]) admits solutions ([7]). The first form 
gives an opportunity to relax the regularity requirements for weak solutions to the MAE, such that the 
integral in the l.h.s. of the identity ([SJ is well-defined. The third form is an integro-differential equation of 
the first order, which might be useful for development of an algorithm for computation of a solution ([7]) 
using transformations of spatial variables. This paper is mostly concerned with the Fourier integral form 
of the MAE used to prove existence of a small-amplitude solution of the form ([7]) for a weakly varying 
r.h.s. in ([T|) and to formulate an algorithm for computation of a solution ([7]) to IT} in an odd-dimensional 
space. In a test application to a three-dimensional MAE with the r.h.s. reminiscent of mass transportation 
problems considered in cosmology, we have demonstrated that a solution to the MAE with a smooth positive 
r.h.s. can be efficiently obtained by two versions of this algorithm, ACPDM (the algorithm with continuation 
in a parameter and discrepancy minimisation) and AICDM (the algorithm with improvement of convexity 
and discrepancy minimisation). While the latter is suitable for computation of a convex solution for an 
everywhere positive r.h.s. /, the former requires only (/) ^ 0. However, modification of ACPDM for the 
case of a zero-mean space-periodic r.h.s. is straightforward. 

There are some analogies between our method and the inexact Newton-Krylov solver with precondition- 
ing, used in f21]. In our algorithms, solving ^7} or ([5^ . where the easily invertible polynomial part of the 
MAE is separated out, can be regarded as preconditioning; we believe that it is optimal. Indeed, its design 
is based on the algebraic nature of the MAE and primarily on the positivity and bounds of the kernels in 
the Fourier integral form, a very special — and so far not reported — property of the MAE. Furthermore, 
our algorithms are of the Krylov type. However, we do not rely on Newton iterations: instead of solving a 
succession of Newton problems, each with a complexity comparable to that of the MAE, we tackle the MAE 
directly. 

In this paper, space periodicity and Fourier decompositions play an important role in two respects: 
First, we rely on the Fourier methods to find the optimal coefficients in the equations (P7)l or 
Note that space periodicity is not required for the derivations — we consider Fourier integrals, and not 
Fourier sums, and the algorithms do not require computation of these integrals. Second, our algorithms are 
realised in the spectral form, because we seek space-periodic solutions. In a periodicity domain the required 
computations of the inverse Laplacian, as well as numerical differentiation become trivial, when spectral 
methods are applied. Recall that Fourier methods have also the advantage of being more accurate: when 
the resolution is increased, Fourier series converge to a solution together with all the derivatives that the 
solution possesses, whilst finite differences do not approximate derivatives beyond their fixed order. We 
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would also use the spectral approach, for instance, to solve numerically Dirichlet or Neumann boundary 
value problems for u' (see (O) in a rectangular box. 

However, the geometry of the domain, where the solution is sought, may prohibit the use of spectral 
methods. What happens when the MAE must be solved for more complicated boundary conditions in re- 
gions whose geometry is not rectangular (parallelepiped-shaped)? Our algorithms do not inherently rely on 
spectral methods and will remain applicable. Since we are solving the MAE in terms of the Laplacian of the 
unknown function, we would have to invert the Laplacian with suitable regular boundary conditions. Finite 
differences can be applied in conjunction with our method to carry out this task and for computation of 
derivatives involved in the function F appearing in the r.h.s. of p2p and (j32p . Key equations for our algo- 
rithms, ((27)) and ([32]), are formulated in the physical space and can be used in conjunction with our method 
for a variety of boundary conditions. Their convergence properties remain of course to be investigated — for 
instance, whether the theoretically predicted values (pS)) will still be the optimal choice for the coefficients 
of the polynomial in the l.h.s. of (j32l) . 

An attractive feature of our algorithm is its simplicity: our Fortran-95 source code realising EAICDM is 
below 18 KB (672 lines long, not including the source for the Fast Fourier Transform). 

Hereafter we list some open questions. Under which conditions does the generalised MAE 

N 

{1-P)J2 amV'^+pdet \\V-%,,,^ + S,j\\ = f/{f) (39) 

m=0 

with the coefficients (I25|) . which is solved for a set of p by ACPDM, possess zero- mean space-periodic 
solutions for all < p < 1? Does its discretisation always have a solution, as long as the generalised 
MAE itself does? Does the solution of (15^1) depend analytically on the parameter p, and hence polynomial 
extrapolation for p = 1, that we employ, is mathematically sensible, or should another asymptotics near 
p = 1 be assumed? What is the optimal choice of the sequence of values of p to be used by ACPDM? Which 
scalar product is optimal for acceleration of convergence of stabilised iterations? Can Chebyshev techniques 
(3^ be used to improve efficiency of the iterative processes? 

Additional questions emerging outside the main topic of this paper — numerical methods for solution 
of the MAE — also cannot be avoided, since any information concerning the structure of solutions to 
Monge-Ampere problems of the cosmological type with a high-contrast r.h.s., formulated in Section 5, can 
be incorporated into specialised solvers in order to improve their performance (like it has proved possible to 
accelerate computations about 4 times just by taking into account in AICDM convexity of solutions to the 
MAE with a positive r.h.s.). The questions are: What is the asymptotics of solutions in the small parameter 
6 determining the width of the localised objects? Figures 4 and 5 show that the space is divided into regions 
of dominant influence of each object. Also, notable are almost axisymmetric structures in the Laplacian of 
the solution around the lines connecting centres of objects, seen on Figure 3. Can these geometric features 
be identified by performing the asymptotic analysis of the problem? How does the contrast number measure 
numerical complexity of the MAE and, in particular, the condition number of the linearisation near the 
solution (for a given spatial discretisation)? 
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